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ABSTRACT 

This survey paper describes recent developments in the area of parametrized variational 
principles (PVPs) and selected applications to finite-element computational mechanics. 
A PVP is a variational principle containing free parameters that have no effect on the 
Euler-Lagrange equations. The theory of single-field PVPs, based on gauge functions 
(also known as null Lagrangians) is a subset of the Inverse Problem of Variational Cal- 
culus that has limited value. On the other hand, multifield PVPs are more interesting 
from theoretical and practical standpoints. Following a tutorial introduction, the pa- 
per describes the recent construction of multifield PVPs in several areas of elasticity 
and electromagnetics. It then discusses three applications to finite-element computa- 
tional mechanics: the derivation of high-performance finite elements, the development 
of element-level error indicators, and the construction of finite element templates. The 
paper concludes with an overview of open research areas. 




1. INTRODUCTION 


This paper is the written version of an invited presentation to the First International 
Seminar in Mechanics (SIMG) held at the Institute de Mecanique de Grenoble on May 1992. 
It reviews the recent development and selected applications of parametrized variational 
principles (PVPs) in mechanics. 

A PVP derives from a functional with free parameters if the Euler-Lagrange equations 
turn out to be independent of those parameters. In this survey a sharp distinction is made 
between single-field PVPs, in which only one primary field is varied, and multifield (mixed) 
PVPs. Single-field PVPs are well known but lack practical importance. They are used 
here only in the tutorial introduction of Section 3. Multifield PVPs are less understood 
but far more interesting from the standpoints of theory and applications. 

The study of multifield PVPs was originally motivated by the desire of finding a variational 
basis for some high-performance finite elements, as outlined in Section 5.6. As sometimes 
happens, this modest goal led to unexpected discoveries, chief among them a general 
parametrization of the functionals of classical elasticity. Subsequently more applications 
to finite elements were revealed. 

The paper is organized as follows. Section 2 is philosophical in nature and speculates on the 
potential impact of this emerging area. Section 3 is a tutorial that presents basic concepts 
and definitions using single-field PVPs. Section 4 introduces multifield PVPs. Sections 5 
through 8 outline the construction of PVPs for several application areas, with emphasis on 
linear elasticity. Those readers with interest restricted to theoretical aspects may skip Sec- 
tions 9 through 12, which present applications of PVPs to finite elements in computational 
mechanics. Of these the first application: development of high-performance elements, is 
the most developed one. The other applications: error estimation and templates, are in 
an exploratory stage and rely largely on conjectures and numerical experimentation. The 
paper concludes with an overview of topics that may deserve further research. 

2. PHILOSOPHICAL DIGRESSION 

Necessity is the mother of invention. Although this aphorism is normally applied to tech- 
nological and industrial advances, there have been instances of its validity in mathematical 
physics. Three examples can be cited. 

Operational Calculus. Invented by Oliver Heaviside to solve the ODEs and PDEs brought 
forth by the initial applications of Maxwell’s electromagnetic theory in the 1880s. Heavi- 
side’s calculus was based on heuristic rules, most of which were eventually justified by the 
development of transform theory and associated axiomatics. By now operational calculus 
is a standard branch of applied and theoretical mathematics. 

Delta Calculus. Although the underlying ideas existed in latent form prior to 1920, delta 
functions were popularized by Paul Dirac as a tool to formulate one of the three versions 
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of the “old” quantum mechanics. Dirac’s recipes were eventually justified twenty years 
later by the theory of distributions and generalized functions. Delta calculus is presently 
an important analytical and computational tool for engineers and scientists. 

Finite Element Method. The pioneering 1943 paper by Courant [15] attracted no attention 
then because no need was addressed and digital computers had not arrived. The rapid 
development of the FEM in the aerospace industry can be traced back to efforts in the 
early 1950s to formulate a satisfactory computer-based procedure for the stiffness analysis 
of sweptback, fixed delta wings [61]. This involved an extension to continua of by-then 
established discrete-element procedures collectively known as “matrix structural analysis” 
[29,34,35,36], which had culminated in the elegant unification of Argyris and Kelsey [3]. 
The method was named in 1960 by Clough [14]. In this early period primary concern 
was given to matching experimented results rather than to the underlying theory. The 
connection with variational principles and approximation theory was established in the 
1960s after the FEM proved successful on a much wider class of problems. 

A common thread runs through these examples: 

1. An urgent application need is established. 

2. Response to need produces an unpolished tool. 

3. Tool is recognized as useful beyond its initial purpose. 

4. With interest from mathematicians aroused, theory (and eventually axiomatization) 
follows. 

As regards parametrized variational principles, only the first two steps have been verified. 
The motivating need was the justification of high-performance finite elements. Although 
there are indications of usefulness beyond that context, such applications remains to be 
firmly established. In any event, much theory is lacking. 

3. PARAMETRIZED VARIATIONAL PRINCIPLES: CONCEPT 

3.1 A Single-Field Example 

Let u = u(x) G C 2 [a,b] be a real- valued function, while jd is a free parameter. Consider 
the parametrized functional 

II(u;/?) = ±/ F{u-p)dx = \! [-(t/) 2 + 2/W + u 2 ] dx, (1) 

J a J a 

where ( )' = d( )/dx. The Euler- Lagrange equation associated with <511 = 0 is 

E{u) = u" + u = 0, (2) 

which is independent of /3. The statement 611 = 0 will be called a parametrized variational 
principle , or PVP. Notice, however, that the natural boundary conditions 

F ul (a) = (3u(a) - u'(a) = 0, F U '(b) = f3u(b) ~ u'(b) = 0, (3) 
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are not independent of /? unless u(a) = u(b ) = 0. 

This example can be easily generalized. Chose ah arbitrary function Q(u , x ) 6 C 2 [a, 6] and 
consider 

n(«; $) = i / [-(„') 2 + 2/3Qu' + u 2 ] dx, (4) 

*/ a 

The Euler- Lagrange equation for (4) is gain (2) — the simpler form (1) is in fact obtained 
if Q — u. Another generalization (noted by a reviewer) is 

/ b 

|_( u <m))2 + -(- (u^- 15 ) 2 ] dx, u (j) = d m u/du m , (5) 

for m > 0, whose Euler- Lagrange equation is u^ m+2 ^ + = 0. Despite their simplicity 

these functionals have engineering significance; e.g. the suitably normalized (5) with m = 2 
governs the buckling of an Euler-Bernoulli beam. 

3.2 Principal and Gauge Functionals 

The explanation for the behavior of the preceding functionals is straightforward. Consider 
for simplicity (2) in which the parametrized term is separated as 

n( w; /?) = np(u) + /?n G ( u ), (6) 

IIp(u) = | f [— (u') 2 + u 2 ] dr, n G (u) = f uu' dx = | [u 2 (6) — u 2 (a)] . (7) 

J a J a 

Here Hp is the principal functional whereas II G is a gauge functional. The latter is also 
called a null Lagrangian in the literature. Because II G depends only on boundary values, 
its contribution to the Euler- Lagrange equation of II is obviously zero thus cancelling the 
effect of (3. 

The converse statement is also a well known theorem of variational calculus: if the Euler- 
Lagrange equation of a single-field functional vanishes identically, that functional must 
depend only on boundary values. See, for example, Section 1.18 of Fox [30]. 

The second example (5) illustrates the addition of an arbitrary one- dimensional gauge func- 
tion. Should the single-field functional depend on several independent variables r,y, ..., 
one can obviously add the divergence or rotor of multidimensional functions multiplied by 
free parameters, because application of Green’s or Stokes’ theorems reduces such functions 
to boundary terms. 

3.3 Terminology 

A functional that contains one or more free parameters, such as n in the foregoing example, 
is called a parametrized functional. If its Euler- Lagrange equation is independent of the 
parameters, the stationarity condition <5n = 0 is called a parametrized variational principle 
or PVP. 
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relations 


Fig. 1. Diagram sketching Strong, Weak and Variational Forms, and 
relationships between form pairs. (Weak Forms are also called 
weighted-residual equations, variational equations, Galerkin 
equations and integral statements in the literature.) 

A PVP is most useful from the standpoint of applications if the value of the functional, 
evaluated at an extremal, is independent of the free parameters. This value has often the 
meaning of energy. Such a principle will be called an invariant parametrized variational 
principle or IPVP. The example functional (3) yields an IPVP if [ u"(b )] 2 = [u"(a)] 2 or 
if [u'(6)] 2 = [u'(a)] 2 . (To prove the latter, insert u = —u" into YIq and integrate.) If 
[u(6)] 2 = [u(a)] 2 , II is independent of /? for any i/(r), not just an extremal. This will be 
called an absolutely invariant PVP, or AIPVP. 

The most useful IPVPs and AIPVPs, however, are those in which invariance does not 
depend on boundary conditions. Such functionals will be encountered in following sections. 

3.4 Connections with the Inverse Problem 

The study of single-field PVP’s constitutes a particular topic of the Inverse Problem of 
Variational Calculus: given a system of differential equations — herein called the Strong 
Form or SF — find the Lagrangians that have that system as Euler- Lagrange equations. 
These Lagrangians (if they exist) collectively embody the Variational Form (VF) of the 
problem. The Weak Form (WF) of the problem, which is also known by the alternative 
names listed in Figure 1, is an intermediary between SF and VF. Relations between SF, 
WF and VF are annotated in Figure 1. 
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The Inverse Problem linking ordinary differential equations to single-field functionals is 
treated in several monographs [55,57,62]. On the other hand, the multifield case is much 
less developed. Aside from the developments presented here, parametrized functionals 
have also been studied by the Beijing school [37,38] but without establishing connections 
to the Finite Element Method. 

4. MULTIFIELD VARIATIONAL PRINCIPLES 


Single-field PVP have limited practical and theoretical value. After all, gauge functions 
have been known for a long time but never attracted much attention. Multifield (mixed) 
PVPs axe more interesting and fruitful, both theoretically and practically. In this section 
two examples are worked out to unveil the flavor of the more complicated cases and expose 
the reader to the notational conventions used in following sections. 

4.1 One-Dimensional Example 


As our first encounter with a multifield PVP, the following 6-coefficient generalization of 
(2) to two independently varied fields: u and p = it', is postulated: 




j 11 

j\2 

j 13 


3 12 

J 22 

j '23 

1 

.713 

]23 

jzz _ 


u 

u' ► 
P , 




z T Jz dx, 


( 8 ) 


in which 



\ 


u 


J 11 

J 12 

3 13 


V j = 

j 12 

J 22 

j 23 

p J 


J 13 

3 '23 

333 _ 


( 9 ) 


are the generalized field vector and the functional generating matrix , respectively. This 
parameter matrix — the “kernel” of the quadratic form in the z vector — may be taken as 
symmetric because only its symmetric part participates in the first variation. This general 
notational arrangement is followed in more complicated linear problems of mathematical 
physics presented later. 


The Euler- Lagrange equations for 6 U II = 0 and ^ p II = 0 are 

Eu : j 11 « + il 2 «' + ilZP - (j 12 U 1 + j 22 u" + j 2 z)p' = 0, 
E p : ii3« + J 23 W' + JzzP — 0. 


( 10 ) 


Consistency of (10) with the field equations u 1 — p = 0 and u + p' =0 dictates that J be 
of the form 

a j3 0 

J= | 0 0 -a\ . (11) 

0 —a a 


Thus the functional is found to depend on two independent free parameters: a and (3. 
Is this PVP invariant? At an extremal, p = u ' , which replaced into (8) yields II = 
^ J^(au 2 — fiu' 2 )dx. Therefore this is not generally an IPVP. 
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A more satisfying parametrization that automatically verifies invariance is obtained by 
taking three independent fields, as in the example that follows. 

4.2 A PDE Example: 2-D Isotropic Poisson’s Equation 


As second example we consider the isotropic Poisson’s equation (the Laplace equation with 
a source term) posed over a finite two-dimensional region ft bounded by a curve T: 


kV 2 u 


/ d 2 u d 2 u\ 

\jW + lh?) 


= ~f 


in ft, 


( 12 ) 


where u = u(x,y ) is the unknown function, k(x,y ) > 0 and f(x,y ) are given scalar 
functions in ft. Boundary T is decomposed into T u UT g , on which Dirichlet and Neumann 
(flux-type) boundary conditions, respectively, are imposed: 


u = u on r u , q = k — = k (grad u) T n = q on T q , (13) 

where u and q are prescribed on T u and T ? , respectively, and n is the exterior unit normal 
on r. The single-field functional II(ti) associated with (12) and (13) is well known: 


n = U(u) - P(u), 




q u dT. 


( 14 ) 

Here U(u ) has the meaning of internal energy whereas P(u ) is an external energy associated 
with the source and prescribed-flux terms. As noted in Section 3.2, functional (14) can be 
trivially parametrized by adding multiples of the divergence or rotor of gauge functions. 
But such PVPs have no practical importance. 


To begin the construction of a multifield PVP we introduce the two intermediate vector 
fields: gradient g and flux vector p, as candidates for independent variation: 

8 _ {fc}“ g " d “"{*/B’ p = { p' } = * g = srad “ = A } ' 

(15) 

The PDE (12) decomposes into the three field equations 


g = Vu = grad it, p = kg, V T p + / = div p + / = 0, (16) 

where div = V T is the divergence operator. In mechanical applications these are called 
the kinematic, constitutive and balance equations, respectively. 

The three field equations (16) and two boundary conditions (13) collectively make up the 
Strong Form (SF) of the isotropic Poisson’s equation. This SF is graphically represented 
in Figure 2 using a modified Tonti diagram. Relations such as g = grad u are called 
strong connections (which means that they are enforced point by point) and depicted as 
solid lines. The main departure of Figure 2 from Tonti’s original diagrams [47,59,60] is the 
explicit separation of field equations and boundary conditions; this has been found useful 
in teaching variational methods. In addition, a graphical distinction is made between 
unknown and data fields, as indicated in Figure 2, also for instructional reasons. 
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Unknown field 


Data field 


Strong connection 


Fig. 2. Graphical representation of the Strong Form (SF) of the 
isotropic two-dimensional Poisson’s equation (12) with 
the B.C.s (13) as a modified Tonti diagram. 


4.3 A 3-Field PVP for Poisson’s Equation 

The configuration of (14) suggests trying to parametrize the internal energy U as a three- 
field quadratic form, leading to the arrangement displayed below in full: 



i \ 

Pz 

T 

'in 

in 

in 

7 12 

713 

in' 


' k~ l Px ' 

Py 


in 

in 

i 12 

j 12 

7' 13 

7' 13 


k~ l P y 

kg x 

> 

i 12 

in 

722 

722 

723 

723 

< 

9z > 

kg y 


j 12 

ii2 

7 22 

722 

723 

723 


9 y 

k du/dx 


ii3 

iis 

7 '23 

723 

733 

733 


du/dx 

k k du/dy , 


J 13 

ii3 

7 23 

723 

733 

733 . 


. du/dy , 


(17) 


in which for simplicity the explicit dependence of U on the j coefficients is dropped from 
its arguments. The symmetry of the kernel matrix can be justified as in the previous 
example, whereas its 2 x 2 block structure is a consequence of avoiding distortions in the 
vector- component contributions to the internal energy. This functional is fully specified if 
the 3x3 generating matrix J, which has the same form of (9), is given. 

Now (17) looks unduly complicated for such a simple problem. At the same time, what 
is being varied is not easily seen. We clarify and simplify this form in two steps: passing 
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to matrix-vector notation, and then applying the primary- versus derived-field convention 
explained below: 


U(u, g,p) = \ f 
Jn 

=ij 

Jn 



ini J13I 
ii 2I J22I J23I 
J13I J23I J33IJ 



- X T 

P 1 

ini ii 2 i ii 3 i 

K1 

< P 9 \ 

ini j2 2I J23I 


P“J 

.J13I J23I J33I. 

lg"J 


(18) 


> da. 


Here I denotes the 2x2 identity matrix, p 9 = kg, p u = k g u = k grad u, g p = k *g, 
g“ = grad u. This notational convention, introduced in [19,20], is based on two rules: 

1. A varied (primary) field is marked with a superposed tilde such as u or p. This allows 
to reserve tildeless symbols such as u or p for generic or exact fields. 

2. A derived (secondary) field is identified by writing its “parent” primary field as su- 
perscript; e.g. p u = k grad u is the flux associated with the varied field u. 

Of course at the exact solution of (12) all p’s and g’s coalesce, but the distinction is crucial 
in variational-based approximation methods. 

Note the pleasing appearance of the last term in (18): the notation groups fluxes on the 
left and gradients on the right. Flux times gradient is internal energy, so the kernel matrix 
simply weights, through the j coefficients, the nine possible combinations p T g p , p T g, . • . 
etc. (It is possible to further streamline (18) into U = I z T Wz da, as in the last of 
(8), but this is too compact for most developments.) These notational conventions are 
especially helpful in the more complicated application problems presented in Sections 5-7. 

The first variation of U may be compactly expressed as 

SU = J [(g) T <5p + (p) T £g - (div p) T <$uj dV + J (p) T n8udS, (19) 

in which g, p and p denote the combinations 

g=jllg p +jl 2 g + jl3g“, P= J 12 P +j 2 2P ff +J23P U , P=il3P +;23P 3 + J33P U - (20) 

On linking 6U with the variation of the parameter-free external- energy term P given in 
(14), we obtain the Euler-Lagrange equations in a 


E p : g= 0, Eg : p= 0, E u : div p +/ = 0, 


( 21 ) 


while the Neumann boundary condition on T 9 is <J= (p) T n = q. Consistency with the 
three field equations (16) and the boundary condition q = q leads to constraint conditions 
on the j coefficients. These can be expeditiously obtained by noting that at the exact 
solution of the Poisson’s problem, p = p = p 5 = p u and g = g = g^ = g“. Consequently 


in + in + ii3 — o, in + ii2 + ii3 — o, in +ii2 + ii3 = 1. 


( 22 ) 
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It follows that the functional ( 18 ) combined with P yields a PVP, which is in fact a three- 
parameter family. If we chose the diagonal entries in, j 22 and j 33 as the free parameters, 
the constraints (22) are met if the functional generating matrix takes the form 



in 2 ( — in ~ J22 +i33 — 1) 

2 (-in +i22 - 

- jzz + 1) 


ill ~ s 3 —$2 

J = 

J22 

§ (in -i22 - 

j 33 + 1) 

= 

~ s 3 ]22 —-Si 


symm 

j 33 



. —$2 —Si j'33 


( 23 ) 


Here the negatives of the three off-diagonal elements of J axe abbreviated to .s 1 , s 2 and S3 
for use below. 

The choice j 33 = 1, others zero, yields the well known single-field functional ( 16 ). Other 
choices for J axe discussed in conjunction with the classical elasticity problem in Section 
5 , which has a similar parametric structure. 

Using the decomposition of J as sum of rank-one matrices 



'0 

0 

o' 


'0 

0 

o' 


1 

0 

-1" 


1 

-1 

o' 

J = 

0 

0 

0 

+ Si 

0 

1 

-1 

+ s 2 

0 

0 

0 

+ S3 

-1 

1 

0 


0 

0 

1_ 


0 

-1 

1_ 


-1 

0 

1 


0 

0 

0_ 


one can rewrite the three Euler-Lagrange equations (21) in a form that illuminates the 
weighted-residual connection to the field equations ( 16 ): 

E P ■ ^ 2 (g p - g u ) + s 3 (g p - g) = 0, 

E g : 5 s(p 9 -p) + si(p fl -p u ) = 0 , ( 25 ) 

E u : div [p“ + s : (p u - p 5 ) + s 2 (p u - p)] + / = 0 . 

What happens if, say, s 2 = s 3 = 0 ? Then E p becomes an identity and g drops out as 
an independently varied field [observe that in =0 because of (22) and consequently the 
first row and column of J vanish.] Similarly if si = s 3 = 0, E g becomes an identity and 
p drops out as varied field. The case sj = s 2 = 0 reduces E u to div p u + / = 0 but 
three-field principles axe still possible because one may select ju — j 22 = —j \ 2 = p, where 
p is arbitrary; setting p — 0 gives back the simplified functional ( 14 ). 

The form ( 25 ) show that the s’s, or their reciprocals, can be interpreted, as weights on 
the field equations. No such flexibility is available with single-field functionals because 
the parameters factor out at the first variation level. This is the key reason behind the 
importance of multifield PVPs. 

Is this PVP invariant? At an extremal the p’s and g’s coalesce. The internal energy 
reduces to f a p T g dV multiplied by ju +ji 2 + . . . + J33, which is one because of (22). Thus 
we have an IPVP. The same property holds for all PVPs presented in Sections 5 - 8 , and is 
crucial in the application to finite element error-estimation discussed in Section 10. 

Does ( 18 ) embed all possible functionals for the internal energy of the Poisson’s equation? 
No! In this formulation u must be a primary field in n, although p and/or g may drop 
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Fig. 3. Graphical representation of the Strong Form (SF) of the 
primal (displacement-based) formulation of classical 
elasticity. Refer to Figure 2 for display conventions. 


out for certain choices of J as discussed above. To make u disappear, the last row (and 
column) of J must have ail zero entries, which contradicts the last of (22). Consequently 
functionals with internal energy of the form ?7(p,g), U( p) or U( g) escape this framework. 
This point is further elaborated in Section 5.5. 

5. CLASSICAL ELASTICITY 

“Classical elasticity” is short for the more precise “compressible linear hyperelastotatics.” 
This is the application area in which the multifield PVPs discussed here originated in 
response to needs from finite element technology. As a result, it is still the best developed 
one in terms of FEM applications. 

The ensuing discussion emphasizes three-dimensional elasticity. Specialization to struc- 
tural models such as beams, plates and shells, as well as the derivation of canonical func- 
tionals for these cases, is covered from a modern viewpoint by several authors, e.g. Reddy 
[53,54] and Hartmann [31]. 

5.1 Governing Equations 

Consider an body of volume V referred to a rectangular Cartesian coordinate system Xi , 
i = 1,2,3. The body is bounded by surface S of external unit normal n = n,-. The surface 
is decomposed into 5 : Sd^St- Displacements d = d{ are prescribed on Sd whereas surface 
tractions t = f,- are prescribed on S t . Body forces f = /, are prescribed in volume V. 

The three unknown internal fields are: displacements u = u,-, strains e = e,j and stresses 
cr = <Jij. The stress traction vector on S is <r n = <r n , = ajiUj (summation convention 
implied). To facilitate the construction of elasticity functionals in matrix form, stresses 
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and strains axe arranged in the usual 6-component vector forms 

Cr T = [&U <722 <733 <712 <723 <731 ] j 

e = [en e22 ^33 2ei2 2 e 23 2e3i ] , 

These fields axe connected by the kinematic, constitutive and balance equations 


e = Du, <t = Ee, D 7 <t + f = 0, (27) 

where 

'dfdx 1 0 0 ' 

0 86 x 2 0 

D - ° ° d,dx 3 (28) 

“ d/dx 2 d/d Xl O’ ^ 

0 djdx z d/dx-i 

_d/dx 3 0 d/dxi _ 

is the 6x3 symmetric-gradient operator, its transpose D the 3x6 tensor-divergence 
operator, and E is the 6x6 stress-strain matrix of elastic moduli arranged in the usual 
manner. The boundaxy conditions are 

u = d on Sdt tr n = t on S t . (29) 


The field equations (27) and boundaxy conditions (29) make up the Strong Form (SF) of 
the classical-elasticity problem. The SF is graphically represented in Figure 3 using again 
a modified Tonti diagram of the primal (displacement-based) formulation. 

5.2 Weak and Variational Forms 

From the Strong Form one can proceed to several Weak Forms (WF) by selectively “weak- 
ening” strong connections. A variation- homogenization and integration process then yields 
a Variational Form (VF) as was sketched in Figure 1. Two specific functionals of classical 
elasticity are worked out as examples before proceeding to the general parametrization. 

Figure 4 depicts the Weak Form (WF) of the Potential Energy (PE) functional for classical 
elastostatics, in which the displacement field u is the only primary field. The WF diagram is 
adapted from the Strong-Form Tonti diagrams of Figures 2-3, with some additional graphic 
conventions explained in Figure 4. Furthermore the derived field notation introduced in 
Section 4.3 appears again. For example, e u = Du is identified by the “tracer superscript” 
u that remind us that these are displacement-derived strains. 

Two strong connections: the balance equations D T cr + f = 0 in V and the flux boundary 
condition: = t on S t , have been weakened by making them satisfied only in an average 

sense: "■ 

f (D r <r u + f)w dV + f «-t)wdS = 0, (30) 

Jv Js , 
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Varied (primary) field 


Derived (secondary) field 


Data field 


— — Strong connection 
«w Weak connection 


Fig. 4. The Weak Form (WF) of the Potential Energy principle of classical elasticity. 
Also known as the Principle of Virtual Work. The function inner-product 
abbreviations (a, b)v d = J y a T bdV and [a,b]$ d = J s a r b are used 
for labeling the weak connections because f is not available for graphics. 


where w are the “virtual displacements” in the parlance of mechanics. To pass to the VF, 
set w = (5u, homogenize variations and integrate obtaining the Potential Energy functional 

IIpE(u) = f f ( * u ) T e u dV - f f T u dV - f t T u dS. (31) 

Jv Jv Js t 


As second example, consider the derivation of the Hellinger-Reissner (HR) stress- 
displacement functional. Figure 5 depicts the appropriate Weak Form (WF). The varied 
fields are displacements u and stresses d\ Three weak links have been introduced: 


/ (e“ -e")6<TdV + f (D T <r + {)6udV + f (fr n - t)6udS = 0 
Jv Jv Js t 

where cr n = VijUj . From this one can pass to the Hellinger-Reissner functional 


(32) 


IWu,<T)= / dr T (e u — ^e ff )dV — f f T udV — f VudS. (33) 

Jv Jv Js t 

There are variants of both the PE and HR functionals in which the displacement boundary 
condition u = d is only weakly enforced. This results in the appearance of an additional 
Sd boundary term displayed in (35) below. 
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Eg 




Fig. 5. The Weak Form (WF) of the Hellinger-Reissner principle of 

classical elasticity. (For notation explanation refer to Figure 4.) 


5.3 Canonical Functionals and Hybridization 

The PE and HR functionals are instances of the so-called canonical functionals. According 
to Oden and Reddy [46,47] there are seven canonical functionals for the primal formulation 
of linear selfadjoint boundary value problems that fit the SF scheme of Figures 2 and 3. 
For classical linear elastostatics these are listed in Table 1, along with the names of the best 
known ones. (Oden and Reddy list eight functionals but one is a trivial combination of 
PE and HR.) This number expands to 14 if one includes the dual formulation of elasticity 
based on potentials, but such generalization — which is discussed in detail by Oden and 
Reddy [47] — will not be considered here. 

Are there more? The number increases if one allows hybrid functionals into consideration. 
All functionals of classical elasticity can be expressed in the form analogous to (14): 

H = U - P, (34) 

where U is the volume term that represents stored strain energy and the potential P 
includes body-force and boundary contributions. The conventional form of P is: 

P c (u,<r,e)= f i T udV + f (a- n ) T (u-d)dS + [ t T \idS. (35) 

Jv J Sd Js t 


w 
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Table 1. The Canonical Functionals of Classical Linear Elastostaiics 

(after Oden and Reddy [47]) 


Varied fields 

Acronym 

Functional name 

u 

PE 

Potential energy 

cr 

CE 

Complementary energy 

e 


No name 

u, cr 

HR 

Hellinger-Reissner 

u, e 

SDR 

Strain-displacement Reissner 

<r, e 


No name 

u, cr, e 

HW 

Hu-Washizu 


where <r n is generally a weighted combination of tr u , <r e and cr as shown in Eq. (43) of 

Section 5.4. The Sd term drops if the displacement boundary condition u = d is strong, 

as in (31) and (33); in such case P c depends only on u. 

Two more forms appear when one considers hybrid functionals for a finite element dis- 
cretization of V. (Such functionals can also be constructed for bodies with internal phys- 
ical interfaces.) In the following Si is the counted- twice union of element interfaces not 
on S, while d and t are boundary- displacement and boundary-traction fields, respectively, 
that are independently varied on S t : 

P d (u,a, e,d) = f f T udV + ( (o- n ) r (u — d) dS + f t T u dS + f <x£(u — d)dS. (36) 
Jv Js d Js, JSi 

P*(u,t)= f f T udV+ f t T (u-d)dS+ f t T u dS+ [ t T u dS. (37) 

Jv Js d Js t JSi 

Use of either of these potentials allows u to be discontinuous across internal interfaces. 
In (36) fields d and u are weakly connected on S',, whereas in (37) fields t and <r n are 
weakly connected on S,. With these additional choices one can have, for example, three 
variants of the Hellinger-Reissner functional: II^ R = Uhr — P c , Lljj R = U R r — P d y and 
= Uhr ~ P *• The last two were called d-generalized cind f-generalized function- 
als, respectively, in [19,20]. The hybrid functionals of Pian and Tong [51,52,63] can be 
precipitated as special cases of these two forms. 

When these hybrid variants are accounted for, the number of primal canonical functionals 
roughly triples (it is not exactly 7x3 = 21 because certain hybridizations are incompatible 
with some functionals of Table 1). But in fact all of them are mere “points” in the space of 
parametrized functionals constructed below. Therefore the correct answer to the previous 
question is that there is an infinite number of functionals from which classical elasticity 
can be derived. 
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1 0 0 
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I SDR 0 -1 i 
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-1 0 1 

1 1 -1 


1 -1 0 
-1 0 1 
0 1 0 


Fig. 6. Representation of the three-parameter PVP for classical 
elasticity in (jn, j'22, .733) space. Generating matrices 
for interesting functionals are shown near the “points.” 


5.4 Heuristic Parametrization 

The strain energy portion U of II = U — P can be written down as follows for the PE, HR 
and HW canonical functionals: 


Vpe( u) = 1 f 

Jv 


0 0 0 
0 0 0 
0 0 1 


U H r{u,&) = \ f 

Jv 


v 


Uhw( u, e) = | 

Jv 


-I 0 I' 
0 0 0 
10 0 

0 -II 
-I I 0 
10 0 


where I is the 6x6 identity matrix. This pattern suggests trying the parametrized form 


U(u,<r y e) = \ f 

Jv 


v a u 


ini J12I J13I 
J12I J22I J23I 
J13I J23I J33I. 
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This is formally similar to (18) except that the <r and e vectors have 6 components, and 
the generalized field vectors 18 components. 

The first variation of U can be compactly written 


su 


•Jvl& 


8(7 + ( <t) t Se — (div er ) 7 <5u 


dV + 


L^ n) 


T SudS, 


(42) 


in which e , a and a denote the weighted combinations 


e=ine p + j l2 e + j 13 e u , <r= j 12 a- + j 22 <r 9 + j 23 tr u , &= ji 3 d- + j 23 cr 3 + j 3Z (r u , (43) 

Combining 8U with 6P yields <511. For example, if P is taken as P c of (35) 


<5n 




e) T 8(7 + {<t) t 8e — r T <5u 


dV + 


f (r n -t) T 8udS+ [ (d — u) T <5 <r n dS , (44) 
Js t Js* 


where r = div <r -f-f are internal equilibrium violations. Consistency arguments again show 
that the coefficients must satisfy the constraints (22) reproduced here for convenience: 


j 11 +j 12 + Jl3 = 0, j 12 +J22 +J23 = Jl3 + J23 + j 3 3 = 1- (45) 

This leaves 6 — 3 = 3 free parameters. As usual the 3x3 matrix J of entries jkt is called 
the functional generating matrix. The special settings 



'0 

0 

o' 


'-1 

0 

r 


' 0 

-1 

r 


'0 

0 

o' 

J = 

0 

0 

0 

, J = 

0 

0 

0 

, J = 

-1 

1 

0 

, j = 

0 

-1 

1 


_0 

0 

1 _ 


-1 

0 

0 


1 

0 

0 


0 

1 

0 


(«) 

yield the PE, HR, HW and SDR canonical functionals, respectively. Other interesting but 
anonymous functionals in ju, j 22 , j 33 space are marked with a o in Figure 6. A glance 
at this figure derails the lofty status accorded the Hu-Washizu (HW) functional in many 
textbooks as the most general three-field functional of elasticity. 

The weighted field-equation representation analogous to (25) is 


E e : ^(e** — eU ) + S 3 (e** — e) — 0, 

E„ : s 3 (a- e - <x) + si(cr e - er u ) = 0, (47) 

E u : div [<r“ + si(cr“ - a e ) + s 2 (<r“ - a)] + f = 0 , 

where scaling factors Si , s 2 and S 3 are defined in (23). The general Weak Form for arbitrary 
parameters is depicted in Figure 7. As can be observed this has become more complex, 
requiring some study to sort out. The “strength” of the weak connections can be measured 
by the value of the scaling factors as discussed in Section 4.3. If factor pairs vanish “box 
merging” occurs as explained therein. 
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A 

l 1 



p 1 

d 


u 


f : 



Weak connections conjugate to 6a 
Weak connections conjugate to 5e 
Weak connections conjugate to 5u 

Fig. 7. Weak Form for the general parametrization (40). Numbers annotated 
near weak connections are the scaling factors in Eqs. (46). 

5.5 Functionals Without Independent Displacements 

There is a subset of elasticity functionals without independently varied displacements; for 
example the well known principle of Complementary Energy (CE). In Felippa and Militello 
[22] it is shown that all functionals with independently varied stresses can be embedded in a 
one-parameter form that includes CE as a special case. The only functional left out of these 
three-parameter and one-parameter framework is the unnamed canonical functional whose 
only independent varied field is strains. Summarizing: all canonical primal functionals of 
linear elasticity can be covered with a 3-parameter family II(u,dr, e), a 1-parameter family 
II(<x,e), plus a 0-parameter family 11(e). This statement can be extended to encompass 
hybrid variants. 

5.6 Historical Note 

The three-parameter form (41) was not the first one constructed. A one-parameter subset 
that connects points PE and HR of Figure 6 was discovered in 1987 and published in 
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Fig. 8. Generalization of the Strong Form of classical linear elastostatics 
(diagramed in Figure 3) to encompass material incompressibility. 


[19,20], This was the byproduct of an effort to variationally justify the Free Formulation 
(FF) of finite elements [12,13]. The effort was actually motivated by criticism of a follow- 
up paper [18] by a reviewer who justly observed that the excellent performance of the 
FF-based plate-bending element presented therein had computational but no theoretical 
basis. 

Subsequent developments were serendipitous. Another one-parameter form, joining 
“points” HR and HW of Figure 6, had been constructed [40,41] to justify the Assumed 
Natural Strain method [7,32,50,58]. This work evolved later into the Assumed Natural De- 
viatoric Strain (ANDES) formulation [22,43]. Comparison of the SFF and ANDES results 
led to the general parametrization (41) published in [22,23]. 

6. ENCOMPASSING INCOMPRESSIBLE BEHAVIOR 

The functionals of linear elastostatics studied in the previous section fail if the material 
is incompressible. To encompass incompressibility one must begin by splitting the stress 
tensor into deviatoric and mean stresses, and the strain tensor into deviatoric and mean 
strains. In Cartesian tensor notation: 

Vij = Sij - pSij, P = trace a, 

' eij = gij — \B8ij, 0 — — e,, = — trace e = — div u. 

Here p denotes the (actual) pressure, 6 the total strain condensation, 8 t] is the Kronecker 
delta, and the summation convention applies over the range 1,2,3. These tensor compo- 
nents are arranged in 6- vector form following the prescription (26) so we can rewrite the 
splitting (48) as 

<7 — s — ph, e = g — #h, h T = [l 1 1 0 0 0]. 
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(49) 





The kinematic relations and balance equations are still e = Du and D <r + f = 0 as 
in (27). The boundary conditions are still (29). As far as the constitutive equations is 
concerned, the assumed split is 


s = Gg, p = kO, 


(50) 


where G is a (generally anisotropic) deviatoric-stress-to-deviatoric-strain matrix invertible 
for all values of compressibility, and k is thrice the bulk modulus I\ . Incompressibility 
occurs in the limit k — ► oo. Note that the model (50) assumes that the material is 
volumetrically isotropic although deviatoric anisotropy is permitted. The Strong Form 
diagram shown in Figure 8 summarizes the governing equations. 

The governing functionals are again II = U — P. The internal energy U has been 
parametrized by Felippa [24] in the form 


U(u,s,p,g,d) = 1 f 

Jv 


\ 

s 

T 

' in I 

712 I 

713 I 

714 h 

7i5h 

iieh' 

s* 


712 I 

722 I 

723 I 

724 b 

725 h 

726h 

s" 

J 


713 I 

723 I 

733 1 

734 h 

735 b 

73eh 

\ 

p 

r 

7i4h r 

724h T 

734h r 

744 

J 45 

746 

P 9 


7 i5h r 

725 h T 

735h r 

J '45 

755 

756 

Ip* J 


.ii 6 h r 

726h T 

736h r 

746 

756 

766 . 





g 

< 

f 


9 p 


9 


,o u > 


dV, 


(51) 

where I is the 6x6 identity matrix, h is defined in (49), and the derived fields are 
s 9 = Gg, s u = G -1 Du, p 9 = k9, p u = —k div u, g a = Gs, g“ = GDu, 9 P = k~ x p 
and 9 U = — div u. In the incompressible limit 9 P vanishes. The analysis of [24] shows 
that consistency of the Euler-Lagrange equations with the field equations requires that the 
following nine constraints be satisfied: 


3 a +;'i2 +j 13 = 0, 
724 + J25 + 726 = 0, 
746 + 756 + 766 = 1? 


714 +715+716 = 0, 
713 +723 +733 = 1) 
744 + 745 + 746 = 0, 


712 + 722 + 723 = 0, 
734 + 735 + 736 = 0, 
745 + 755 + 756 = 0. 


(52) 


Furthermore, to include the incompressible limit, j 55 = j $$ = 0. Consequently there are 
21 — 11 = 10 free parameters. The simplest generating matrix that satisfies all constraints 


is 


J = 


"in 

i 12 

i 13 

714 

7 15 

716 


'0 

0 

0 

0 

0 

O' 

i 12 

722 

723 

724 

725 

726 


0 

0 

0 

0 

0 

0 

ii3 

723 

733 

734 

735 

736 


0 

0 

1 

0 

0 

0 

ii4 

724 

ju 

744 

i45 

746 


0 

0 

0 

-1 

0 

1 

ii5 

725 

735 

745 

755 

756 


0 

0 

0 

0 

0 

0 

,ii6 

726 

736 

746 

756 

766 . 


.0 

0 

0 

1 

0 

0. 


(53) 


This leads to a generalization of the PE functional (31) that exactly accomodates incom- 
pressibility. Other choices are examined in [24], where stress-strain splittings more general 
than (48) are considered. 
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Fig. 9. 


Strong Form diagram for micropolar elastostatics without couple stresses. 


7. MICROPOLAR ELASTOSTATICS 

The extension to a linearly elastic micropolar medium without couple stresses has been 
investigated by Felippa [25]. In this elasticity model the stress tensor rand the strain tensor 
7 are no longer symmetric. They are decomposed as r = <r + s, and 7 = e -J- <f>, where cr 
and e are the symmetric parts and s and <f> the antisymmetric parts, respectively. (Symbol 
s does not have the meaning of Section 6.) General, symmetric and antisymmetric tensors 
that arise in this model may be arranged as 9, 6, and 3-component vectors, respectively, 
to facilitate the use of matrix notation in the parametrized functionals. Notational details 
are given in that reference. 

The rotational portion of u now enters the field equations. The infinitesimal rotation 
vector (also called macrorotations) is u> = Ru, where R is the rotor or curl operator. 
The microrotation vector 6 is introduced so that the antisymmetric strains are given by 
<f> = uj — 0. The irrotational part of u is u e , which generates the symmetric strains 
e = Du = Du e , where D is again the symmetric gradient operator. Finally, body- volume 
actions are decomposed into f = b + c, where b and c are the body force and body couple 
per unit volume, respectively. Classical elasticity is recovered if c vanishes. 

The balance equations Eire D^tr + f = 0 and R^s + c = 0. The constitutive equations 
are <x = Ee and s = G <j>, which may be merged as r = C 7 , where C is block diagonal. 
These field equations, plus the boundary conditions u = d on Sj and r n = t on St, are 
summarized in the Strong Form diagram of Figure 9. 

The functional II again decomposes into U — P. The parametrized form of the internal 
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energy U has the structure 


U(u, er, s, e, 



r <x ' 

T 

7 'iile 

712I6 

713I6 

<r e 


712I6 

722I6 

723I6 

<r u 

► 

713I6 

723I6 

733I6 

s 

0 

0 

0 



0 

0 

0 



0 

0 

0 


0 0 O' 


' e * 7 ' 

0 0 0 


e 

0 0 0 


e" 

744I3 745I3 746I3 

< 


745I3 755I3 756I3 


4 > 

746I3 756I3 766I3 . 


± u9 

4 > 1 


( 54 ) 


where 16 and I 3 denote the identity matrices of order 6 and 3, respectively. The derived 
quantities axe cr e = Ee, cr u = EDu, s ^ = G<£, s u0 = G(Ru — 0 ), <p s = G -1 s and 
(j> ue = Ru - 0 . 


The zero entries of the kernel matrix in (54) reflect the orthogonality of symmetric and 
antisymmetric parts. The analysis carried out in [25] shows that consistency with the field 
equations requires that 


in +ii2 + in — 0, ii2 +J22 +J23 — 0, iia + 723 +733 = 1, 

735 + 745 + 746 = 0? 745 + 755 + 756 = 0> 746 + 756 + 766 = !• 


(55) 


Consequently there axe 12 — 6 = 6 free parameters in (54). One possible choice for the 
generating matrix that meets these constraints is 
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1 
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0 
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to 
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. 0 

0 

0 

1 

0 

0. 


(56) 


An alternative application of these functionals is to non-polar media in which c = 0 
although the symmetry of the stress tensor is only enforced weakly. Such functionals 
furnish a basis for construction of finite elements for classical elasticity with independently 
varied rotation fields, a subject that has recently been given impetus by the paper of 
Hughes and Brezzi [33]. In fact the choice (56) yields one of the canonical functionals 
derived in that paper. The Appendix of [25] extends these PVPs to micropolar media with 
couple stresses. 


8 . ELECTROMAGNETODYNAMICS 

As final example, consider classical electromagnetodynamics based on Maxwell’s equations. 
In this model the unknown internal fields are the magnetic and electric fluxes H, and D 
the magnetic and electric intensities B and E, the magnetic potential A and the electric 
potential $. The prescribed volume fields are the current intensity j and the electric source 
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Fig. 10. The Strong Form of isotropic linear electromagnetodynamics, with 
boundary conditions omitted for simplicity. Here e and p denote 
the electric permittivity and magnetic permeability, respectively, of 
the medium. Other symbols are defined in the text. A superposed 
dot denotes derivative with respect to time t. 


p. The field equations axe summarized in the SF depicted in Figure 10, in which boundary 
conditions have been omitted for simplicity. 

The general functional has the form II = U — S — B, where U is the generalized elec- 
tromagnetic energy, S includes source terms and B contains boundary-closure terms. A 
parametrized form of U studied by Felippa and Schuler [28] is 
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} dVdt, (57) 


where I is the 3x3 identity matrix and c is the speed of propagation of EM disturbances. 
(The parametrization coefficients for this functional are denoted by g instead of j to avoid 
confusion with the symbol for current density.) The derived terms in (57) are H B = B, 
H a = p- 1 curl A, D E = eE, etc. The analysis of [28] shows that the g coefficients must 
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satisfy the 12 consistency constraints: 


011 + 0 12 + 013 

= 0, 

014 + 015 + 016 = 0, 

012 + 022 + 023 

= 0, 

024 + 025 + 026 = 0, 

013 + 023 + 033 

= 1, 

016 + 026 + 036 = 0, 

034 + 035 + 036 

= 0, 

046 + 056 + 066 = ~ 1, 

014 + 024 + 034 

= 0, 

044 + 045 + 046 = 0, 

015 + 025 + 035 

= 0, 

045 + 055 + 056 = 0> 


( 58 ) 


which leaves 21 — 12 = 9 free parameters. The simplest choice for the generating matrix is 
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which, on including prescribed sources p and j, yields the two-field Lagrangian 


(59) 



1 1 V x A] 


-i e 


v$ + 


dA 

~dt 


-3 T A-p§ \ dVdt. 


(60) 


Felippa and Schuler [28] extend this PVP to the situation in which the current density j 
is unknown, in which case Ohm’s law has to be adjoined to the governing equations. This 
extension is important as a departure point for the finite element analysis of Type I and 
II superconductors. 


9. APPLICATIONS TO FINITE ELEMENT METHODS 

From previous examples it is seen that multifield PVPs are interesting in their own right 
because they provide a continuous space of functionals. Canonical functionals are “instance 
points” of this space. Moving in this space is equivalent to changing weights on the field 
equations. A theorem proven for a parametrized functional is “economical” in the sense 
that it need not be redone for specific instances. 

Traditionally books dealing with applications of variational calculus in mechanics jump 
from one specific canonical functional to another. This approach has two instructional 
disadvantages. First the reader is left with the impression that only a finite number of 
functionals exist. Second, what ought to be shared properties are proved over and over, a 
boring repetitive feat that can easily span hundreds of pages. 

Aside from this theoretical and educational unifying value, multifield PVPs offer applica- 
tions to variationally-based methods of approximation in general and the Finite Element 
Method (FEM) in particular. Following are three intriguing subjects pertaining to FEM, 
which in the author’s opinion merit further exploration. 
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(A) Improved Approximation. Approximate FEM solutions generally depend on the 
free parameters while the converged solution (assuming a convergent approximation 
method) does not. Can the parameters be chosen to improve the quality of approxi- 
mation on a fixed finite element mesh? 

(E) A-posteriori error estimation. Assume the PVP is an IP VP. Suppose we obtain FE 
solutions for two different parameter sets on the same mesh. As the mesh is refined the 
values of the functional for both solutions approaches each other over each element. 
Can this property be exploited to construct local (element level) error estimators? 

(T) Templates. Using PVPs, can we derive universal “templates” for FE matrix expres- 
sions that covers all convergent finite elements of a given freedom configuration? 

The answers to these questions are in different state of development. 

Exploiting (A) leads to high-performance finite elements. This is an application that has 
been investigated over the past four years with focus in classical elasticity, and has led 
to production-level finite elements. Knowledge as to the effect of parameter selection is 
increasing, but much remains to be done. 

Exploiting (E) leads to element-level error estimators. Only limited numerical experimen- 
tation has been carried out, and most of the theory remains to be developed. 

The idea of templates (T) has been explored for some specific finite elements, but it remains 
largely a conjecture. 

The next three Sections describe the state of these applications in further detail. 

10. APPLICATION 1: HIGH PERFORMANCE FINITE ELEMENTS 

A high performance finite element (HPFE for short) for structural mechanics is a simple 
element that delivers results of engineering accuracy on coarse meshes. 

By “simple” it is meant that the element has few degrees of freedom, all physical, preferably 
at corners only. Engineering accuracy typically means less than 1% error in displacements 
and 5% in stresses. The term “coarse mesh” is subjective. In terms of structural applica- 
tions, it characterizes a mesh that an experienced modeler would use to capture the basic 
physics with a minimum number of degrees of freedom. This number is of course problem 
dependent: a coarse mesh for a simple supported, uniform loaded plate would typically be 
a2x2or4x4 mesh with less than 100 freedoms; on the other hand a coarse mesh for a 
complete aerospace vehicle may involve 50,000 to a million freedoms. 

Mandatory attributes for a HPFE are: convergent, frame invariant and rank sufficient. 
Desirable attributes include: 

• Yields similar accuracy in displacements and stresses. 

• Is relatively insensible to geometric distortion. 

• Is mixable with other elements. 

• Provides effective a-posteriori error estimators to drive mesh adaptation procedures. 
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• Fits naturally into displacement-based FE programs. 

• Extends readily to nonlinear and dynamic problems. 

The goal of attaining reasonable accuracy with coarse meshes of simple elements is the key 
one, however, and deserves further comment. This requirement should not be confused with 
fast asymptotic convergence for fine meshes. Simple elements cannot effectively compete 
with higher order elements in this regard, and are not effective in applications that demand 
very high accuracy. What matters is the accuracy obtained for the meshes typically used 
in complex engineering systems, as discussed above. 

10.1 Unification 

Many tools and techniques have been developed over the past three decades to develop 
HPFEs. The most practically important are: (a) incompatible shape functions; (b) hybrid 
and mixed formulations; (c) reduced, selective and directional numerical integration; (d) 
assumed natural strains and the Free Formulation. These four techniques originated in 
the mid 1960s, late 1960s, early 1970s and early 1980s, respectively. 

The approach recently pursued by the author has relied upon the use of PVPs with hybrid 
treatment of element interfaces. More specifically, three ingredients are used: 

1. The parametrized internal energy functional (41). 

2. A d-generalized hybrid treatment of interfaces through the forcing potential (36). 

3. Additional assumptions (sometimes called “variational crimes” or “tricks”) that can 
be placed on a variational setting through Lagrange multipliers. 

The first two ingredients says that HPFEs can be effectively constructed using hybrid 
PVPs of the form 

n d (u, cr, e, d) = £7(u,cr, e) — P d (u,a-, e,d), (61) 

where U and P d are given by (41) and (36), respectively. Observe that P d contains indepen- 
dently varied displacements d on element interfaces. This allows the internal displacements 
u to be nonconforming and in some cases (ANDES elements) ignored entirely. 

Throughout this section we shall deal with individual elements only. Element identifiers, 
which will be required in Section 11, are omitted for simplicity. 

For mechanical elements modeled by classical elasticity, the HPFE derivation method 
summarized in Figure 11 leads to finite element stiffness equations that decompose as 

(K b + K a )v = p, (62) 

where v are visible degrees of freedom (those matched with other elements, also known 
as connectors ), p are associated node forces, Ki is the basic stiffness matrix, which is 
constructed for convergence and mixability, and K* is the higher order stiffness matrix, 
which is constructed for stability and accuracy. This decomposition was first obtained, 
without variational arguments, by Bergan and Nygard in the derivation of the unsealed 
Free Formulation [12]. A key result of the multifield variational formulation is that only 
K/, depends on the functional parameters [22,23]. 
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Fig. 11. Summary of assumptions made in constructing a high-performance mechanical 
element from the PVP (61). For the basic stiffness an interior constant stress 
field o- = & and an interelement-compatible boundary displacement field dj are 
assumed. For the higher order stiffness either an internal displacement field u 
is assumed with the FF and its variants, or an internal strain field e is assumed 
with the ANS formulation and its variants. (Only the deviatoric part of e 
appears in the ANDES variant.) The boundary displacement field d>, need 
not be the same as d&, although it often is. In the FF and variants thereof, 
d/i and u are collocated at node points. 


10.2 The Parametrized 3-Node Bar Element 

The 3 -node bar element depicted in Figure 12 will be used to illustrate the steps in the 
derivation of HPFE stiffness equations from PVPs. The element has total length 2 h 
and constant axial rigidity EA. The applied force per unit length is /( x). The internal 
fields are the axial displacement u(z), axial strain e(x) = u'(x), and internal axial force 
N(x ) = Aa(x ) = EAe(x). The 3 degrees of freedom (dof) axe the node displacements iq, 
V2 and v 3 along the bar axis. For an internal finite element of this type, the d-generalized 
parametrized functional ( 61 ) reduces to 



where N e = EAe, N u = EAu ' , e N = (EA) 1 N, e“ = u 1 , and iV= —s 2 N — siN e +733^“. 

To construct the element we make assumptions on the varied fields N, u, e and d. In 
expressing u and e it is convenient to introduce the hierarchical displacement of node 2: 
v 2 = v 2 — 4. V3 ) ? which is the midpoint deviation from linearity. We also introduce 

the natural coordinates = |(1 — x/h) and £3 = |(1 + x/h) for convenience. 
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Fig. 12. Three-node bar element used as example of derivation of 
stiffness equations from the ^-generalized principle (63). 

The internal force is constant over the element: N = N. The mean force N is a degree of 
freedom that (thanks to the hybrid treatment) can be eliminated at the element level. 

The axial displacement is expanded into a basic (linearly- varying) component which 
represents rigid body modes and constant strain states, and a higher-order (quadratically- 
varying) component u 

u(x) = U b (x) + Ufc(x), U b (x) = UjCi + v 3 ( 3 , u k (x ) = 4^2 C 1 C 3 - (64) 

The assumed strain is also split into a basic component (the average or mean strain e) and 
a higher-order component called the deviatoric strain: 

e = e + e h , e = (v 3 - Vi )/h, e h (x) = y[//sign(x) - 2(1 - (65) 

Observe that e b is represented here by a constant-plus-linear function that is odd in x, and 
contains a coefficient [i that weights the “mixing” between the constant and linear parts. 
If n = 0 the strain varies linearly as in the quadratic- displacement element and e = e u . 
If fi = 1 the strain is constant over subelements 1-2 and 2-3 and agrees with that of the 
“macroelement” constructed with two linear 2-node bar elements. The assumption (65) 
satisfies the compatibility conditions f° h e b (x)dx = — f 0 e b (x)dx = V 2 for any //. 

Finally, the boundary displacements are identified with the end-node displacements: 

d(—h ) = d = —v\, d(h) = d - v 3 . (66) 

— h h 

Because u(x) also agrees exactly with uj and v 3 at the end nodes, the boundary term in 
(63) vanishes. Inserting the above assumptions into (63) produces a quadratic form in 
the degrees of freedom v\, V 2 , v 3 and N. Rendering the form stationary with respect to 
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these freedoms leads to a system of linear equations that immediately gives N = EAe = 
EA(v 3 — v i ) / h. Elimination of N yields the element stiffness equations 



0 -1 

0 0 

0 1 


+ 0 


2 EA 


0 0 
1 0 
0 0 


ui 

V2 



(67) 


where forces p \ , p 2 and p 3 result from the nodal lumping of f(x) and 


0 = |[( 4 -m)Jh + (^ 2 - M)j22 + W33 + (4 - /i)] = \[p 2 s\ + 4s 2 + (p 2 — 2p + 4)s 3 + 4], (68) 
Expressing (67) in terms of the total midpoint displacement v 2 rather than v 2 yields 


Kv = (K b + K h )v = 


l EA 

' 1 

0 

-1' 

E A 

l 2 h 

0 

0 

0 

+ 0— 

-1 

0 

1 

2h 


EA 

6h 


3 + 3 0 -60 -3 + 30 

-60 120 -60 
-3 + 30 -60 3 + 30 



(69) 


Here we see the emergence of the stiffness decomposition (62). The basic stiffness Kj 
contains no parameters and has rank one, so by itself it is rank deficient. The higher order 
stiffness has also rank one and its addition stabilizes K as long as 0 > 0. 

Although the FEM solution converges for any positive 0 , maximum accuracy in smooth 
problems is obtained for 0 — 4/3. This yields the well known 3-node quadratic bar 
element. [To check this, take j 33 = 1, j u = j 22 = 0 in (68) to get the PE functional, which 
annihilates the separate strain assumption (65).] On the other hand, setting 0—1 yields 
the 3-node “macroelement” equivalent to assembling two linear bar elements. [To check 
this, take j 22 = 1, ju = / 33 = 0 in (68) to get the HW functional and set p. = 1]. 

The foregoing derivation appears to be overkill for such a simple element. Indeed (69) can 
be obtained without recourse to any variational principle. Direct physical reasoning (as in 
the original FF), algebraic arguments relying on subspace decomposition of node motions, 
or even finite differences: all methods lead to (69). Such universality is the motivation 
behind the template conjecture discussed in Section 12. In two and three dimensions and 
especially complicated elements, however, variationally-based methods do have a place. 

10.3 The Parametrized 2-Node Plane Beam Element 


As a second one-dimensional example consider a two-node element for an Euler-Bernouilli 
(C 1 ) plane beam of constant rigidity El and span h. The four nodal degrees of freedom 
are the transverse end displacements w j and w 2 and the end rotations Q\ and 6 2 as shown 
in Figure 14. Any paxametrization method, whether variationally based or not, produces 
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Fig. 13. Two-node Euler-Bernouilli beam element. 


a decomposed stiffness matrix of the form 


K = K b + K h = 


El 


'0 

0 

0 

0 ' 


4 

—2h 

-4 

-2 h' 

0 

1 

0 

-1 

a EI 

+ p h 3 

—2h 

h 2 

2h 

h 2 

0 

0 

0 

0 

-4 

2h 

4 

2h 

,0 

-1 

0 

1 _ 


-2 h 

h 2 

2h 

h 2 , 


(70) 


where /? > 0. Both Kj and K& have rank one, and combine to provide the proper rank of 
two. The (optimal) Hermitian-cubic beam element is recovered if (3 = 3. 

An interesting property is that the Timoshenko ( C ° ) beam has exactly the same stiffness 
decomposition but the optimal (3 is then 3 + 3 6(EI/GA)h~ 2 , where GA is the section- 
averaged shear rigidity; if GA —* oo or h — ► 0, /? — ► 3. The distinction made in the copious 
finite element literature on C° versus C l beam elements is seen to be artificial. 


10.4 Three High-Performance Triangles 

For multidimensional elements the general parametrized functional (61) has not been yet 
exploited in its full generality. As of this writing only two one-parameter subsets have 
been studied in some detail: 

1. A stress-displacement d-generalized functional n^(u,cr, d) associated with the Free 
Formulation (FF) and two variants thereof: the Scaled Free Formulation (SFF) in- 
troduced in [13] and the Extended Free Formulation (EFF) proposed in [21]. Its free 
parameter is called 7. Reduces to PE for 7 = 0 and to HR for 7 = 1 (see Figure 14). 

2. A stress-strain-displacement d-generalized functional n£(u, <r, e,d) associated with 
the Assumed Natural Deviatoric Strain (ANDES) formulation. Its free parameter is 
called a. Reduces to HR for a = 0 and HW for a = 1 (see Figure 14). 

The most successful multidimensional HP elements constructed to date using these subsets 
are depicted in Figure 15 and briefly described below. The EFFAND element is a plane 
stress triangle (membrane component in shell element) with nine degrees of freedom, three 
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Fig. 14 . The one-parameter functionals associated with the Free Formulation 

(FF) and Assumed Natural Deviatoric Strain (ANDES) formulation of 
finite elements, depicted in the general (ju, j’22, .733) space of Figure 6 . 


of which are corners “drilling” rotations 6 Z . It represents the culmination of a development 
began by Bergan and Felippa [13] with the Seeded Free Formulation. The latest element was 
derived with both the Extended Free Formulation [21] and the Assumed Natural Deviatoric 
Strain (ANDES) formulation [2,26,27]. Triangles optimized for bending behavior were 
found to coalesce, and this optimal element named EFFAND. 

The insensitivity of EFFAND with respect to element aspect ratio is illustrated in Table 
2, which pertains to the beam-under-pure-bending test problem shown in Figure 16. The 
results of Table 2 are taken from Felippa and Alexander [27]. In this Table ALL-xx denote 
three numerically integrated versions of the 1987 Allman rank-sufficient triangle [1], CST 
is the constant strain triangle [61], and FF is the first-generation triangle of this type 
constructed by Bergan and Felippa [13]. 

AQR is a Kirchhoff (thin) plate bending element constructed by Militello and Felippa 
[43]. The element has the standard nine degrees of freedom and is based on the ANDES 
formulation. An extensive set of numerical tests shows that AQR generally outperforms 
other elements of this type tested to date as regards distortion sensitivity. 

Finally, the EFFAND-AQR shell element is a combination of EFFAND (membrane com- 
ponent) and AQR (bending component) to form a triangular-facet shell element. It has 18 
degrees of freedom: 9 for membrane plus 9 for bending. This element has been found to 
excel in modeling complete aerospace vehicles for stiffness and vibration analysis because: 

(a) Membrane/bending accuracy is balanced. 
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Fig. 15. Three high-performance elements constructed with PVPs: the 
EFFAND membrane 9-dof element, the AQR plate-bending 
9-dof element, and the EFFAND-AQR 18-dof shell element. 


(b) Global results — e.g. frequencies — of engineering accuracy may be obtained with 
“coarse” meshes (for a complete aircraft, 50,000 to 200,000 dof). This is important 
for full-vehicle aeroelastic simulations. 

(c) The presence of a complete set of rotational freedoms eliminates juncture modeling 
problems. 

(d) The element simplicity helps implementation on medium and fine-grained parallel 
supercomputers with limited memory per processor. 

11. APPPLICATION 2: LOCAL ERROR ESTIMATION 

Is 'i J'»! fci. 3 • • ' - ■ 

As second PVP application we take up the subject of a posteriori estimation of the local 
discretization error of a computed finite element solution. These estimators may be used 
to guide computer-driven mesh adaptation processes based on r, h or p techniques. 

The theory of FE local estimation and mesh adaptivity has received substantial attention 
over the past 17 years since the appearance of the original papers by Babuska [4] and 
Sewell [56]. Much of the theory pertains to conforming finite element models based on the 
Potential Energy principle or its equivalents. Because elements based on these assumptions 
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Fig. 16. Test problem for hi gh -performance membrane element: 16:1 

cantilever of rectangular cross section under pure end moment M . 
Poisson’s ratio v = 0, E adjusted so that exact tip deflection is 
100. A 32 X 2 mesh is shown. Each rectangular macroelement is 
made up of four “crisscrossed” membrane triangles. (After [27].) 


Table 2. Tip Deflections (exact=100) for Beam of Figure 16 


Element 

V 

Mesh: ^-subdivisions x 
32 x 2 16 x 2 8x2 

y-subdi visions 
4x2 2x2 

ALL-3i 

0 

87.99 

75.47 

37.01 

5.51 

0.42 

ALL-3m 

0 

81.02 

51.62 

9.64 

0.74 

0.04 

ALL-7i 

0 

85.43 

67.44 

23.65 

2.55 

0.17 

CST 

0 

53.33 

33.33 

13.33 

3.92 

1.02 

EFFAND 

0 

100.00 

100.00 

100.00 

100.00 

100.00 

FF 

0 

100.25 

99.15 

98.38 

98.08 

97.98 


are rarely the best performing ones, a curious situation arises in which the well developed 
error theory applies to rather uninteresting elements. 

11.1 Error Equations for Poisson’s Problem 

We review the basic error equations for a conforming finite element discretization of the 
Poisson’s problem discussed in Section 3.2 based on the PE-like principle (14). The analysis 
follows essentially Babuska and Miller [5] as well as Babuska, Duran and Rodriguez [6]. 

The bounded domain ft with boundary T is replaced by a regular FE mesh with nonover- 
lapping elements numbered e = l,...N e . Element e has interior ft c and boundary T 6 . 
The latter is composed of N ‘ element sides identified as j = 1, . . . N*. For simplicity 
the geometric discretization error is neglected so that ft = U e ft e and T = U ej T^ V Tj € T. 
Further it is assumed that when an element side belongs to T, it is completely contained 
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in either T u or r ? . The union of element boundaries T® not on T is called r,, which thus 
collects all internal sides. An internal side that is shared by two adjacent elements e and 
d is called r e< *. 

Let u = u(x,y) be the exact solution of the Poisson problem (11)— (12) with / £ L 2 {Tt), 
q € L 2 (r ? ) and u = 0 on IV Let u = u(x,y) £ H 1 ( Q) : = 0 be a conforming finite 

element solution, and e u = u — u the approximation error. The localization of u over 
element e is u e . The source residual is r = V 2 u -f / in each Q e and is conventionally zero 
over element sides. For each element e we denote the boundary-normal flux derived from 
u as q ue = kdu e /dn e where n e is the unit outward normal on T e . On each internal side 
T ed shared by elements e and d we select a common normal direction n ed and define the 
flux jump as 

lq ued ] = q™ cos(n e , n de ) - q ud cos(n d , n de ), (71) 

a value that is independent of the choice for n ed . Let w = w(x,y) £ H J (f2) : ixj = 0 be 
a test function. Introduce the usual bilinear form associated with the internal energy (14) 

U(u,w) = f (div u) T div w dQ,. (72) 

Ju 

On integrating by parts one finds that the error satisfies 


U( 


u — u,w) = U(e u ,w) = J fw dQ + J V 2 u w dQ, + J q w dT — J 

= f rw dfl + f (q — q u )w dT — j [g“ e<i ]uMfr, 

Jn Jr, , Jv^GTi 


du 

dn e 


w dT 


(73) 

which may be called the Weak Form of the Error (WFE). It displays three contributions. 
The first term shows the effect of the source-residual error r = V 2 u + / over each element. 
The second term brings the effect of violating the flux boundary conditions. The third 
term gives the contribution of interelement flux jumps. 


The error energy measure is obtained by taking w = u — u = e“ : 


R d = U(e u ,e u ) = f re u dQ + / (q - q u ) e u dT - Y] [ lq ued ]e u dT, (74) 

Jn Jr , e d J r^er. 

which may be called the Error Energy Equation (EEE). Relations of this nature are the 
basis for developing a posteriori local error estimators for the Laplace and Poisson equa- 
tions as well as (with suitable extensions) the Stokes and elasticity problems. There is an 
abundant and rapidly growing literature in this topic. The state of the art is typified by 
the research of Oden and coworkers [16,48,49]. 

Several technical difficulties, however, may be noted. 

1. Both WFE and EEE involve the exact solution u , which is of course unknown. Con- 
sequently a higher order estimation is inevitable. One way to achieve this is by 
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Fig. 17. Bar problem used to illustrate behavior of HOE error indicator. 


considering element patches (macroelements) and smoothing a, as in the widely used 
projection method of Zienkiewicz and Zhu [67]. Another group of techniques relies on 
hierarchical approximations [64,65,66]. 

2. Many interesting problems display jumps ( e.g . at shocks, wavefronts, plastic bands) 

in the exact solution. The resulting lack of smoothness — in the previous case, u can 
be no longer assumed in — requires substantial revisions in the analysis. 

3. The FE solution has been assumed to be variationally consistent with Eq. (14), that 
is, fully conforming. But in many practical problems, especially in structural and 
continuum mechanics applications, interelement continuity is lost. This may be due 
to several factors. First, as noted in Section 10 most high-performance finite elements 
(HPFEs) abandon conformity from the outset. Second, the treatment of complicated 
3-D structures such as shells, folded plates or stiffened panels often introduces discon- 
tinuities even if the constituent elements (by miracle) were fully conforming. 

4. Although a WFE such as (73) can be formally extended without difficulty to PVPs 
thus fitting HPFEs, the resulting EEE is generally indefinite and thus unsuitable as 
an error measure. It is preferable to keep error measures positive definite, as in the 
example that follows, by concentrating on a single field of interest to the user. 

11.2 The 3-Node Bar Element Revisited 

To motivate the approach discussed below, it is useful to consider again the 3-node bar 
element of Figure 11, now as the problem depicted in Figure 17. This example connects to 
the Poisson’s equation discussed above by simply taking k = EA, u as axial displacement, 
and / as prescribed axial force. Domains O and T reduce to x-aligned segments and end 
node points, respectively. The bar ends are assumed fixed: v\ = v$ = 0; therefore i> 2 = v 2 . 
End fixity is not an important restriction because any element must exactly capture the 
linear motion uj of (64). 

We investigate here an assumed strain element variant of that considered in Section 10.2. 
It is based on the ANDES functional (already mentioned in Section 10.4) that is obtained 
by taking jn = a — 1, j 22 = js 3 = 0, and which reduces to HW for a = 1 and to HR 
for a = 0. The assumptions on N, e and d are as before. The displacement field is again 
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u = Ub + Uh but its higher order part is now 


Uh = v 2 


^sign(i) -(1 -p)|j + 1 


(75) 


rather than (64). This is obtained by integrating the assumed strain (65). 

The applied force is /(x) = /, a constant. The finite element matrix equations reduce to 
the scalar equation k 22 v 2 = p 2 , in which 


P2 = ^ g — hf, k 22 = (3= ±a(4 - 2(* + i?), v 2 = 

The error energy equation is derived from the computed strains: 


4 — jx h 2 f 

2o(4 — 2/j. + fj. 2 ) EA 

(76) 


R = h f 

J-h 


EA 


e(x) — e(x) 


1 2 


dx 


8(1 — ar) 2 (2 — //) + (1 — 2c* + Aa 2 )/j. 2 h? f 2 


12a 2 (4 — 2 fj. + /i 2 ) 


EA 


(77) 


[Using e u yields the same result here because e = e“ on account of (75). Note that only if 
a = 1 and /i = 0 the element yields the exact solution for this problem.] Now comparing 
(77) to the higher order energy (HOE) absorbed by the element 


U h = iv r K,v = \k 22 v\ = 


22^2 = \? 2 V 2 ~ 


(4 - /i) 2 /l 3 / 2 

a(4 — 2fj. + fi 2 ) EA 


(78) 


we note the same dependence on problem properties EA, h and /. The error coefficients 
of R and Uh match if we choose 


a — 
V = 


48 - 24 /i + 3 n 2 ± (4 - /x) v 'W--40/i^"7/2 2 
8(4 — 2 ft + fj. 2 ) 

4-12 a + 4a 2 ± 4\/3o\/3a - a 2 - 1 


(given /i), 


1 — 3o + 4a 2 


(given a). 


(79) 


A particularly interesting combination is a = 1, /i = 2\/3 — 2. If one departs from these 
pairs, Uh and R do not display the same error coefficients. In the accepted FEM 
terminology, (78) would be an error indicator but not an estimator. Of course, should one 
select a and // beforehand, matching is always possible by choosing 3 differently for the 
higher order stiffness and error energy computations. This matching game can be played 
with other functionals and displacement- strain assumptions. 

The example is instructive and may be extended to more general force distributions, ge- 
ometry and boundary condition data, as well as to C 1 and C° beams. But of course 
one-dimensional problems are of little interest in practice. Can the underlying idea of 
“higher error energy” be useful as an indicator or even estimator for more general situa- 
tions? The PVP framework provides the tools to at least pose this problem. 
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11.3 Local Error Estimation Using IPVPs 

Consider again the IPVP of classical elasticity with in which three free parameters are 
generically shown a s a, 0 and 7 in the argument list: 

n(u, or, e; a, (3, 7 ) = U(u, or, e; a, 0, 7 ) - P(u, <x, e), (80) 

Here P stands for either P c , P d or P*. If not P c an additional boundary field d or t would 
appear in the argument list for P. This is of no consequence to the argument that follows 
because the potential P is not affected by the parametrization. 

Now suppose that (80) is used to construct a finite element discretization, while keep- 
ing one or more parameters free. Practically this means that such parameters are kept 
as arguments of the dement stiffness formation subroutines. The discrete approximants 
u(a,/?> 7 )i a(a,(3, 7 ) and e(a,/3, 7 ) obtained on a fixed mesh can be viewed as functions of 
the free parameters. 

Because (80) is a IPVP, U takes the same value if the varied fields u, a, e were set to 
the exact solution fields u, < 7 , e, independently of a, and 7 . Hence, we may expect that 
the difference between two values of U obtained for two different sets of parameters may 
provide a measure on how far we are from the converged solution. That is, the difference 
may be adopted as an error indicator. The most interesting feature of such an indicator is 
that it naturally provides an element level measure, as the following argument shows. 

To tackle the general case first, suppose that the two sets of parameters are (a*, di» 
71 ) and (ao, /?o, 7 o)- The corresponding approximate values for displacements, strains 
and stresses obtained with a given mesh are (ui, 07 , ej) and (uo, d'o, eo), respectively. 
Denote by U{ = U c (u 1 ,d- 1 ,e 1 ,ai,/ 3 i, 7 i) and Uq — U e (u 0 , o'o, eo, Po, 7 o) the value 
of the generalized strain energy evaluated over the e th element of that mesh. Then the 
element error indicator is defined as the difference 

R e = Wt — Uq\. (81) 

The definition (81) apparently requires that the problem be solved twice for a given mesh. 
Additional assumptions and manipulations discussed by Beltran [ 8 ], Militello [42] and 
Militello and Felippa [45] show that the error indicator (81) can be expressed in the one- 
solution form 

U’h = i(V) r KJv*, (82) 

where v e denotes the element node displacement vector extracted from the complete dis- 
crete solution for a single parameter set (a, f3, 7 ). Because only one solution is involved, 
we shall refer to this as a one-solution indicator. Physically (82) is the higher order energy 
(HOE) absorbed by the element and already encountered in Section 11 . 2 . 

Inasmuch as for the HP elements discussed in Section 10 the matrices of the decomposition 
(62) are usually available in separate form, this indicator can be readily computed in an 
element-by-element manner. No calculation of interelement stress-flux jumps is necessary. 

It should be observed that the indicator (82) is related to that heuristically proposed by 
Melosh and Marcal [39] in the context of the SED (Strain Energy Density) method. 
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11.4 Historical Note and Discussion 

Certainly (82) is computationally attractive, as in principle it eliminates all of the tech- 
nical objections raised in Section 11.1 against the “residual plus jumps” estimators. It 
is also physically transparent because it “filters out” contributions from the rigid body 
and constant stress states. Thus, in problems whose analytical solution consist of uniform 
stress states the indicator vanishes over each element, in accordance with the fact that 
any mesh of complete elements should solve those cases exactly. This statement may be 
extended to problems solved by constant stress states that jump at internal interfaces, if 
such interfaces are modeled correctly by the FE mesh. 

The key question is: can £7j£ be considered an error estimator? Little theory is available to 
answer this question, aside from rather trivial one- dimensional cases such as that discussed 
in Section 11.2. Evidence to date is based largely on numerical experiments. 

The first experiments were conducted by Beltran [8] and presented by Beltran and Alarcon 
[9,10,11]. They exploited the Scaled Free Formulation functional II 7 to derive a HOE 
error indicator for plane stress and plate bending problems. This indicator was applied 
to h (element subdivision) and p (element enrichment) mesh adaptation processes. It was 
concluded that performance was generally satisfactory although connection with the actual 
error was erratic. 

Militello [42] conducted experiments in plate bending and shell problems with elements 
formulated with the Assumed Natural Deviatoric Strain (ANDES) formulation. The in- 
dicator was applied to both r (node relocation) and rh mesh adaptation processes. The 
general conclusion is that the indicator was effective in driving those processes in flat 
plates. In shells the performance was less consistent possibly because (82) cannot account 
for geometric discretization error without information from adjacent elements. 

Most experiments so far have dealt with fairly simple elements in which the higher-order 
stiffness is an incomplete-polynomial correction to the basic behavior. This is in agreement 
with the philosophy that HP elements should be “simple but not simplex.” If the HOE 
indicator is to be extended to more refined elements, the appropriate decomposition of the 
element stiffness matrix is not K e = Kj + but 

K‘ = K; + K‘ kc + KJn, (83) 

where the complete higher order stiffness K£ c is produced by displacement-stress-strain 
expansions complete up to a certain polynomial order, and the hierarchical higher order 
stiffness is due to additional terms. The HOE indicator is then ^(v e ) T K. e hh v e . Such refined 
elements appear in p mesh-adaptation processes. 

The connection of such indicators to hierarchical error estimators is an open research topic. 
From evidence gathered to date it appears that free parameters for HOE evaluation need 
to be adjusted from those used for the stiffness calculations — as in the example of Section 
11.2 — so as to get an error estimator rather than just an indicator. In any case, it is 
obvious that if K hh vanishes the HOE indicator is vacuous. Thus complete elements such 
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as the Linear Strain Triangle or Quadratic Strain Triangle [17] should be viewed as useless 
for this purpose. 

12. APPLICATION 3: FINITE ELEMENT TEMPLATES 

The following stiffness matrix decomposition has been observed to recur in several FF and 
ANDES elements (in particular, those mentioned in Section 10): 

K = K 6 + K* = cLEL t +j3TlST h , (84) 

where c is the inverse of the element volume, area or length, j3 is a positive scaling factor, 
L and T& axe geometric matrices, E is the matrix of elastic moduli, and S depends on the 
geometry and constitutive behavior as well as variational principle parametrization. (The 
element index e has again been omitted from (84) for simplicity.) 

Additional properties axe as follows. The geometric matrix L, which is called the force- 
lumping matrix in accordance with the FF terminology [12], is responsible for satisfaction 
of the patch test as well as that of element mixability properties; a detailed discussion of 
these points has been recently given by Militello and Felippa [44]. The geometric matrix Tj, 
annihilates the rigid-body modes and constant strain states. Together with S it provides 
the correct rank to K&, and thus is responsible for stability and accuracy. 

Both L and S may contain free parameters. Any free parameters in L, however, must 
be the same for all elements in an assembly; otherwise convergence and mixability are 
impaired. On the other hand, parameters in S, as well as j3, may vary from element to 
element as long as the correct rank is preserved. 

The form (84) will be called a finite element template, or template for short. The template 
conjecture asserts that (84) is the most general form for all convergent FE of given degree - 
of- freedom configuration. 

A parametrized expression that provides all such elements is called an universal template. A 
parametrized expression that provides a practically important subset of all such elements is 
called a generic template. Equations (69) and (70) provide examples of universal templates 
for the 3-node bar and 2-node plane beam element, respectively. 

Verification of this conjecture for more general elements would have the following projected 
implications: 

1. The confusion associated with the enormous number of “named” elements claimed in 
the literature clears up: such elements, just like the canonical functionals discussed 
previously, become points in a parametric continuum. Among other benefits, such 
unification should help to streamline teaching of finite element theory and hence release 
instruction time to cover problem modeling aspects. 

2. Construction and testing of parametrized FE libraries can be greatly simplified. 

3. Another realm of mesh adaptivity opens up: c adaptation, in which programs can 
pick up elements from the parametrized continuum to strive for maximum accuracy 
on a fixed mesh. 
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4. Several persistent mysteries may be cleared once and for all: convergence, mixability, 
accuracy, locking, spurious modes, and distortion sensitivity. These properties could 
be traced to components of the template (84) and interpreted physically. 

13. CONCLUSION 

The paper has surveyed the development of Parametrized Variational Principles and se- 
lected applications to the Finite Element Method. It is shown that multifield PVPs provide 
a framework by which classical variational principles may be embedded into a parametriz- 
able continuum of functionals. 

Three FEM applications have been reviewed: construction of high performance finite el- 
ements, local error estimation, and element templates. These are in different stages of 
development. While the application to HPFEs has succeeded beyond initial expectations, 
the development of element-level error estimation has primarily relied on numerical exper- 
iments, and the notion of templates is largely a conjecture suggested by isolated examples. 

In the author’s view, the following areas of variational methods merit further development: 

• Parametrization of boundary terms, and use of such functionals in interfacing finite 
element and variational Boundary Element discretizations (or, more generally, Trefftz 
methods). 

• Extension to dynamic and nonlinear problems. It is well known that finite elements 
optimized for linear static analysis may not perform as well when time-dependency 
and nonlinear effects are included. Functional parametrization opens up possibilities 
of extending the aforementioned FEM applications to such cases. 

The following FEM-application areas would benefit from focused research work: 

• Establishment of a general convergence theory of FEM based on PVPs. Areas that 
require special attention are: dependence of element performance on choice of param- 
eters, mixability properties, and limitation principles. 

• Development of a sound a-posteriori error estimation methodology for two- and three- 
dimensional nonconforming elements; in particular establishing connections between 
HOE and hierarchical error estimators. 

• Verification or disproval of the template conjecture. If proven, proceed to exploit the 
benefits outlined in Section 12. 
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